
# Load packages

  library(tidyverse)
  library(haven)
  library(stargazer)
  library(broom)

# Prepare data
  
  # Load dataset
  dt <- read_sav("ObsStudy_00_data_NZ.sav") %>%
    filter(Finished==1, quota_out!="1") 
  
  # Vote 
  vote_labour <- as.numeric(dt$ptv_1)
  
  # Coalition evaluation
  Z <- dt %>% select(contains("coal_rat")) 
  rat_coal <- Z %>% select(-coal_rat_2) %>% 
    mutate_all(as.numeric) %>%
    mutate_all(funs(scales::rescale(.,to = c(0, 1), from = c(1,7)))) %>%
    as.matrix()
  

  # Coalition likelihood
  coal_lik <- dt %>% select(contains("coal_lik_labour"))  %>% 
    mutate_all(as.numeric) %>%
    mutate(sum_exp = exp(coal_lik_labour_1) + exp(coal_lik_labour_2) + exp(coal_lik_labour_3) + exp(coal_lik_labour_4),
           coal_lik_labour_1 = exp(coal_lik_labour_1)/sum_exp,
           coal_lik_labour_2 = exp(coal_lik_labour_2)/sum_exp,
           coal_lik_labour_3 = exp(coal_lik_labour_3)/sum_exp,
           coal_lik_labour_4 = exp(coal_lik_labour_4)/sum_exp) %>%
    select(-sum_exp) %>%
    as.matrix()
  
  # Mean-Variance Model
  N <- nrow(coal_lik)
  
  EV <- function(gamma,Z){
    gamma %*% Z
  }
  
  Vari <- function(gamma,Z){
    gamma %*% ((Z - as.numeric(EV(gamma,Z)))^2)
  }
  
  M <- sapply(1:N, function(i) EV(coal_lik[i,], rat_coal[i,]))
  V <- sapply(1:N, function(i) Vari(coal_lik[i,], rat_coal[i,]))

  # Data Frame
  d_nz <- data.frame(
             "ptv" = vote_labour, 
             "rat"=as.numeric(dt$party_rat_2),
             "sex" = dt$gender,
             "age" = dt$age,
             "edu" = dt$educ,
             "lottery_mean" = M,
             "lottery_variance" = V
             )
  d_nz$ptv <- scales::rescale(d_nz$ptv, to = c(0, 1), from = c(1,7))
  d_nz$lottery_variance <- scales::rescale(d_nz$lottery_variance, to = c(0, 1), from = c(0,0.25))
  
# Saving model results
  
  m1_nz_tidy <- tidy(m1_nz <- lm(ptv ~ lottery_mean + rat + sex + as.factor(age) + as.factor(edu), d_nz))
  m2_nz_tidy <- tidy(m2_nz <- lm(ptv ~ lottery_mean + lottery_variance + rat + sex + as.factor(age) + as.factor(edu), d_nz))

  NZ_2020_OS_Labour_Variance_Estimate <- m2_nz_tidy %>% filter(term=="lottery_variance") %>% select("estimate")
  NZ_2020_OS_Labour_Variance_SE <- m2_nz_tidy %>% filter(term=="lottery_variance") %>% select("std.error")
  NZ_2020_OS_Labour_Mean_Estimate <- m2_nz_tidy %>% filter(term=="lottery_mean") %>% select("estimate")
  NZ_2020_OS_Labour_Mean_SE <- m2_nz_tidy %>% filter(term=="lottery_mean") %>% select("std.error")
  